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Abstract 

We propose a working hypothesis supported by numerical simulations that brain networks evolve based on the 
principle of the maximization of their internal information flow capacity. We find that synchronous behavior and 
capacity of information flow of the evolved networks reproduce well the same behaviors observed in the brain 
dynamical networks of Caenorhabditis elegans and humans, networks of Hindmarsh-Rose neurons with graphs given 
by these brain networks. We make a strong case to verify our hypothesis by showing that the neural networks with 
the closest graph distance to the brain networks of Caenorhabditis elegans and humans are the Hindmarsh-Rose 
neural networks evolved with coupling strengths that maximize information flow capacity. Surprisingly, we find 
that global neural synchronization levels decrease during brain evolution, reflecting on an underlying global no 
Hebbian-like evolution process, which is driven by no Hebbian-like learning behaviors for some of the clusters 
during evolution, and Hebbian-like learning rules for clusters where neurons increase their synchronization. 


Author Summary 

The study of the function of the brain is of primordial importance in neuroscience. Several brain models have been 
studied so far that take into account higher level functions of the external stimulus and the behavioral response 
attributed to ensembles of neurons of cortical areas. If the brain learns by maximizing the Mutual Information 
between stimuli and response, or by updating the internal model of probabilities by using Bayesian techniques, 
or by minimizing the free-energy, it provides little insight about the dynamical mechanisms appearing in brain 
networks when evolved based on such principles. In our work we propose a working hypothesis supported by 
numerical simulations that brain dynamical networks evolve based on the principle of the maximization of their 
internal information flow capacity, i.e. the upper bound for the information transferred per time unit between any 
two nodes. We make a strong case to verify our hypothesis by showing that the neural networks with the closest 
spectral graph distance to the brain networks of Caenorhabditis elegans and humans are the Hindmarsh-Rose 
neural networks evolved with coupling strengths that maximize information flow capacity. We also find that 
synchronous behavior and capacity of information flow of the evolved neural networks reproduce well the same 
behaviors observed in the brain dynamical networks of Caenorhabditis elegans and humans. Finally, we find that 
global neural synchronization levels decrease during brain network evolution. 


Introduction 

A plethora of phenomena in nature can be effectively described by networks. Neuroscientists have used tools for 
the analysis of complex networks that help realize even more deeply the functionality and structure of the brain. It 
was found that many aspects of brain network structures are typical of a wide range of non-neural or non-biological 
complex networks [l ,|2]. One of the main findings in neuroscience is the modular organization of the brain, which 
in turn implies an inherent parallel nature of brain computations [l]. Modular processors have to be sufficiently 
isolated and dynamically differentiated to achieve independent computations, but also globally connected to be 
integrated in coherent functions [I]. It has been revealed that the cortical network is a hierarchical and clustered 
network with a complex connectivity 3 . A possible network description for this modular organization is that 
brain networks may be small-world structured [4| with properties similar to many other complex networks [5j. 
This viewpoint has been driven by the systematic finding of small-world topology in a wide range of human brain 
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networks derived from structural [I], functional 6 , and diffusion tensor MRI [7 studies. Small-world topology has 
also been identified at the cellular-network scale in functional cortical neural circuits in mammals [8] and also in 
the nervous system of the nematode Caenorhabditis elegans ( C.elegans ) [9 . Moreover, this topology seems to 
be relevant for the brain function because it is affected by diseases 10], normal ageing, and by pharmacological 
blockade of dopamine neurotransmission 
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Synchronization is ubiquitous in nature. Insightful findings regarding synchronization in complex networks 
were reviewed recently in Ref. 12 . Recently, synchronization in complex modular or clustered networks has been 
investigated [13 14 . It appears as the interplay between the intrinsic dynamics associated to the nodes of the 
network and its graph topology and connecting functions. In this work, synchronization will be considered as a 
mean to quantify functional behaviors of the brain dynamical networks (BDNs) studied. By BDN we mean a 
network that represents the connectome equipped with neural dynamics for its nodes to account for their time 
evolution. 

Mathematical and computational approaches have a long tradition traced back to the early mathematical 
theories of perception 15 and of current integration by a neural cell membrane [16] . Hebb’s idea on assembly 
formation 17 inspired simulations on the largest computers available at that time (1956) |l8j to understand the 
relation between neural connectivity strength and response, i.e. behavior. It is a learning rule which proposes an 
explanation for the adaptation of neurons during the learning process. The simultaneous activation of pairs of 
neurons leads to pronounced increases in their synaptic strength. There is also the possibility of many other kinds 
of learning 19 . In this work we find evidence of Hebbian-like and no Hebbian-like learning rules characterized by 


the relationship between the addition of synapses and the increase or decay in the synchronization behavior of 
neurons in the evolved BDNs. 

Several brain models have been studied so far that do not take into account particular behaviors of isolated 
neurons, but higher level functions of the external stimulus and the behavioral response attributed to ensembles of 
neurons of cortical areas. There are two classes: Those based on collective functional dynamics of local groups 
of neurons, such as the Wilson-Cowan model for the cortical and thalamic nervous tissue [20], and those based 
on the conditional probabilities (as well as information and mutual information (MI)) of stimuli and responses, 
prominent examples of which are the Bayesian brain hypothesis 21 and the infomax theory 22 . Recently, it was 
shown that most of the probabilistic brain models can be unified under a single free-energy principle 23 , the one 


that interprets the brain as a system that tries to minimize surprises of the sensations from the world. 
If the brain learns by maximizing the MI between stimuli and response 


model of probabilities by using Bayesian techniques 21 , or by minimizing the free-energy 23 
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or by updating the internal 
the surprise of the 


response provides little insight about the dynamical mechanisms appearing in a brain network when it evolves 
based on such principles. The drawback however of the probabilistic brain models is that they describe little about 
the underlying dynamical structure of what is really happening in the neural level and the representation of the 
stimuli 
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For deterministic systems with correlations, an appropriate quantity for measuring the transfer of information 
is the Mutual Information Rate (MIR), the MI per time unit. In Ref. 25 , the authors have developed alternative 
methods to overcome problems that stem from the definition of probabilities from time series and derived an 
upper bound for the MIR, I Cl between two nodes of a complex dynamical network from time averages that do not 
rely on probabilities, but instead on the two largest Lyapunov exponents l \. I 2 of the subspace of the network 
formed by the two nodes (see Eq. (10) in Materials and Methods). In our study, I c stands for the upper bound for 
the information transferred per time unit between any two nodes of the BDN, what represents the information 
flow capacity of the BDN. We discuss more on these in Materials and Methods, Section Upper Bound for MIR. 

Inspired by the infomax theory and by theoretical studies that proposes that the maximization of information 
transmission between subsystems can be used as a principle for understanding the development and evolution 
of complex brain networks (see Ref. [2 and references therein, and Ref. [26]) and, aiming at elucidating the 
microscopic dynamic mechanisms associated to the evolution of brain networks, we propose a working hypothesis 
and, provide evidence that a brain dynamical network may evolve based on the maximization of the information 
flow capacity it can internally handle at each step of its evolution process, i.e. a new inter-neuron connection is 
established if it leads to a subsequent increase of the information flow capacity of the new brain circuitry. Our 
hypothesis is based on the internal neural network dynamics and on a plausible model for brain structure per se 
without the need to resort to probabilistic models based on the external influence in the brain and its response. 

We have been able to show that our evolved BDNs present similar synchronization and information flow 
capacity behaviors with those found for the simulated dynamical networks for the brain structure of the C. elegans 
and humans. Moreover, we show that BDNs evolved with coupling strengths that maximize the information flow 
capacity are the ones with the smallest spectral graph distance from the BDNs of the C.elegans and humans, and 
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that, during the growing process, their information flow capacity increase is related to moderately low amounts of 
global neural synchronization. This work provides ample evidence that brain networks may grow by maximizing 
the capacity of information flow they can internally handle, driven by global no Hebbian-like evolution processes, 
according to which the addition of interconnections between clusters during the evolution process leads to a 
decrease in the global synchronization level of the BDN. This behavior is accompanied by a similar no Hebbian-like 
learning process for neurons in some of the clusters and by Hebbian-like processes for neurons in the remaining 
clusters, leading to an increase in the synchronization level between these neurons during brain network evolution. 
Effectively, the no Hebbian-like mechanism is akin to the unlearning anti-Hebbian mechanism of Crick and 
Mitchison 27 that refers to the elimination of unnecessary connections to prevent overload and to render the 
network more efficient. In our evolution we do not delete links. However, both mechanisms lead to a decrease of 
synchronization and to more efficient networks, being in our case the evolved networks able to maximize their 
information flow capacity. Global synchronization takes into consideration the synchronous behavior of all neurons 
in the BDN whereas local synchronization of neurons in a cluster of the BDNs. Finally, our work shows that 
optimizing information flow capacity leads to evolved networks that are heterogeneous, in accordance with the 
line of research in Ref. [2], where the authors report on a mathematical model for the evolution of heterogeneous 
modules in the brain based on the maximization of bidirectional information flow transmission. 


Results 

In this work we are interested in understanding the relation between synchronization and information flow capacity 
in BDNs constructed by connecting electrically and chemically, Hindmarsh-Rose (HR) neurons (see Materials 
and Methods, Subsection Hindmarsh-Rose Neural Model for Brain Dynamics). Their topologies are given by the 
C.elegans and human brain networks, and those that are the result of an evolutionary process that maximizes 
I c based on interconnected small-world communities. Our intention is not to strictly model the C.elegans and 
human brain neural dynamics but merely to use the topology of their connectomes to study the corresponding 
dynamic brain networks and to compare between the results obtained from these cases. In the following, global 
synchronization measurement (denoted by p) quantifies the synchronous behavior of all neurons of the network 
whereas local synchronization (denoted by p Ci ) quantifies the amount of synchronization between ensembles of 
neurons forming the i-th cluster within the BDN. Both are quantified by the order parameter defined in Materials 
and Methods, Subsection Synchronization Measures in BDNs. 

C.elegans BDN 

We present in panels (A), (B) of Fig. [I] the global synchronization measure p and upper bound of information flow 
I c for the BDN of the C.elegans in the parameter space of chemical coupling g n in [0,2] and electrical coupling gi 
in [0,2]. We refer the reader to Materials and Methods, Subsection C.elegans Date[| A direct comparison between 
the two panels reveals a number of conclusions for the different parameter space regions. At first, for relatively 
high chemical and electrical couplings almost full global synchronization can be achieved (yellow and red regions 
in Fig. [IjA)). For the same region, panel (B) shows an almost absence of capability of information transmission as 
the upper bound for MIR, I c « 0 (dark blue region). High levels of global synchronization accompanied by low 
values of I c indicate that not only neural activities are similar but also they have very low entropy since both 
Lyapunov exponents Ai, A 2 are practically zero leading to their difference I c ~ 0 (see Materials and Methods, 
Subsection Upper Bound for MIR). Secondly, for chemical couplings smaller than 0.3 (i.e. g n £ [0,0.3]) and 
electrical couplings gi £ [0,2], we observe a multitude of different functional behaviors: There are regions of 
high synchronization (red region in panel (A)) and low I c (blue region in panel (B)) and others with exactly the 
opposite behavior (i.e. low global synchronization accompanied by high or information flow capacity). These 
different functional behaviors will become even more evident in Fig. [2] where we plot the regions between the left 
vertical axes and the white dotted lines of Fig. [I] in a finer resolution. 

Our findings suggest that for the C.elegans BDN, for most of the chemical and electrical couplings, global 
brain synchronization is roughly speaking inversely related to the information flow capacity, I c . As we shall see in 
Materials and Methods, Section Similarities between C.elegans and human BDNs for a more detailed analysis of 
the C.elegans and human BDNs, high synchronization as depicted by p implies small information flow capacity. 
The relation between p and I c can be better understood in the basis that typically p oc IA 2 , where A 2 is the second 
largest Lyapunov exponent of the BDN, for the so-called non-excitatory networks [28j|29]. The non-excitatory 
character is expected to be prominent when the electrical coupling has a dominant contribution to the behavior of 
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the network with respect to the chemical, a situation that promotes global neural synchronization. 



Figure 1: Results for the global synchronization and information flow capacity properties for the 
C.elegans and averaged human BDNs. Parameter space for the global synchronization p in panel (A) and 
for the upper bound for MIR, I c , in panel (B) for the C.elegans BDN. Panels (C) and (D) are similar but for the 
averaged global synchronization (p) 6 and averaged upper bound for MIR, (I c )e for the six human BDNs. Here, g n 
is the chemical and gi the electrical coupling of Eqs. §• The regions between the left vertical axes and white 
dotted lines are replotted in Fig. [2] in a finer resolution. 


Human BDNs 

We have performed a similar study for the global synchronization p and upper bound of information flow I c for 
the human brain networks in Materials and Methods, Subsection Human Subjects Data, shown in Fig. 0C) 
and (D). Particularly, we first prepared parameter spaces for all six human subject BDNs with the same coupling 
ranges with those of the first two panels of the C.elegans and then computed their average, presented in Fig. 0 C ) 
and (D). The averaged p and I c quantities from the six subjects are indicated by ( p)q and (I c ) q, respectively. 

A direct comparison between panels (C) and (D) of Fig. [T] for the humans reveals a number of parameter space 
regions associated to different functional behaviors, similar to those observed for the C.elegans. For relatively 
high chemical and electrical couplings almost full global synchronization is achieved (yellow and red regions in 
Fig. [lJC)) since (p)e « 1. For the same region, Fig. [l](D) shows an almost absence of information flow capacity 
as (J c ) 6 ~ 0, the reason being the same as for the C.elegans BDN. For chemical couplings smaller than 0.2 (i.e. 
g n £ [0,0.2]) and electrical couplings gi £ [0,2] we observe regions of high synchronization (red regions in Fig. 
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[lJC) ) and low (J c )6 (blue region in Fig. [IjD)), as well as others with exactly the opposite property, having low 
synchronization and large (I c ) 6- 



Figure 2: Magnification for the global synchronization and information flow capacity properties of 

Fig. 0 Panel (A): Parameter space for p and panel (B) for I c for the C.elegans BDN. Panel (C): Similarly for 
(p )6 and panel (D) for the averaged upper bound for MIR, ( I c )q , of the six human BDNs. Here, g n is the chemical 
and gi the electrical coupling of Eqs. ([2]). 

In order to look deeper into the details of the functional behaviors of these BDNs and, to understand the 
structural and functional similarities between them, we study in Materials and Methods, Section Similarities 
between C.elegans and human BDNs, zoom-in plots of the previous parameter spaces. 


Similarities between C.elegans and human BDNs 

There has been enormous research devoted on the C.elegans worm which has revealed its ability to learn about 
mechano, chemo and thermosensory inputs and stimuli [30,3l]. It was also shown that its neural system has the 
ability to distinguish between tastes, odours or any indication related to the presence or absence of food. It also 
shows different kinds of learning behavior, such as associative (classical conditioning and differential classical 
conditioning), and non-associative forms of learning, such as habituation and dishabituation 32 . These properties 
are reminiscent of the human brain ability to adapt to different stimuli and environments. 

In Fig. [2] we present a finer resolution version of the parameter spaces of Fig. [I] for the C.elegans and humans. 
They allow us to reveal the extraordinary similarity on the functional level, the global synchronization and I c 
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patterns between the C.elegans and human BDNs. In these plots, wherever we observe high synchronization as 
evidenced by orange and red regions in panel (A) for the C.elegans and (C) for the human brain network, the 
upper bound for MIR, 7 C , that stands for the internal information flow capacity of the BDN, is small (blue region) 
and vice-versa. 

Since the C.elegans brain connectivity network is about four times smaller in size than the human BDNs 133] 
studied here, we should expect that the global synchronization and upper bound for MIR patterns could occur for 
different ranges of chemical and electrical couplings as compared to those of Fig. [l] (see for example Ref. 34 ). 
Therefore, a rescaling of the coupling strengths was employed to allow for both BDNs to have the possibility 
to produce equivalent dynamical behaviors. This rescaling is described in Materials and Methods, Subsection 
Rescaling of Chemical and Electrical Couplings for Parameter Spaces of Networks with Different Eigenvalue 
Spectra. It is worth noting that there is an optimal coupling range for both BDNs that allows for large information 
flow capacity in the brain networks (orange and yellow regions in panels (B), (D) of Fig. [5]), a coupling range that 
promotes moderately low global synchronization! 


Functional Properties of the Model for Brain Network Evolution 

In Materials and Methods, Subsection A Model for Brain Network Evolution Based on the Maximization of 
Information Flow Capacity, we propose an artificial brain network evolution model that presents important 
structural and functional properties of the BDNs of the C.elegans and humans. It is based on the combined effect 
of chemical and electrical synapses and, on a topology reminiscent of interconnected brain communities found 
in these BDNs. We use chemical synapses for the communication of neurons of different clusters (inter-cluster 
connections) and electrical for the communication of neurons within each cluster (intra-cluster connections). 

Here, we compare the functional properties of this brain network evolution model with those of the C.elegans 
and humans. By functional we mean the properties of the dynamics of the BDN such as local, global synchronization 
and information flow capacity. Particularly, we report on the similarities we found for the functional properties of 
the BDNs of the C.elegans and humans, and for those of the proposed model for brain network evolution. 

To support further the validity of the presented results for the model for brain network evolution and to show 
its independence on the particular initial small-world cluster configuration, we computed its parameter spaces 
averaging over the functional measurements obtained from five different initial small-world cluster configurations 
(for a discussion about the creation of small-world networks or clusters see Materials and Methods, Subsection 
Analysis of Networks and Communities). In all cases, we used the same evolution process as described in Materials 
and Methods, Subsection A Model for Brain Network Evolution Based on the Maximization of Information 
Flow Capacity. The results of this study are shown in Fig. [3] where (p) 5 stands for the average of the global 
synchronization of the evolved BDNs for the five realizations and (mMIR)s for the average maximal value of the 
MIR for the same realizations. We first provide evidence in panel (A) and (C) to (H) that the finally evolved 
BDN of small-world clusters captures similar local and global synchronization properties to those of the C.elegans 
BDN. The local synchronization p Ci of the i-th community (Fig. [3](C) - (H)) also reproduces similarly the global 
synchronization patterns of Figs. [2](A) and (C) for the C.elegans and human BDNs. Comparing these panels, 
we conclude that the averaged synchronization measure (p)$ of the brain network evolution model (Fig. |3](A)) 
attains almost similar values in the same coupling regions to those of the C.elegans in Fig. HA). The different 
ranges on the horizontal axes of the chemical coupling strength g n can be explained by the fact that the finally 
evolved network consists of 60 neurons whereas the C.elegans of 277 neurons (for more details see Materials and 
Methods, Subsection Rescaling of Chemical and Electrical Couplings for Parameter Spaces of Networks with 
Different Eigenvalue Spectra). 

It is remarkable that the finally evolved BDN exhibits almost identical synchronization and information flow 
capacity features as the BDN of the C.elegans and humans for almost all coupling ranges considered. Particularly, 
focusing on panels (A), (B) (for the C.elegans ) and (C), (D) (for the humans) of Fig. [5]and, on (A), (B) of Fig. [3] 
(for the brain network evolution model), we observe an almost identical pattern of global synchronization and 
information flow capacity as depicted by I c and its averages. Again, here we have used the rescaling in Materials 
and Methods, Subsection Rescaling of Chemical and Electrical Couplings for Parameter Spaces of Networks with 
Different Eigenvalue Spectra, to create networks that can potentially reproduce similar dynamical behaviors for 
the finally evolved BDNs, in agreement with those identified for the C.elegans and human BDNs earlier. 

The upper bound for MIR, 7 C , depends on various factors, such as network topology and connectivity patterns, 
coupling strengths and types (chemical, electrical), synchronicity, etc. The values of 7 C , the information flow 
capacity, can be comparable (or different) for networks with different topologies. This is due to the fact that 7 C is 
a function of the chemical and electrical coupling strengths. The surprising fact is that as we evolve a network 
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Figure 3: Results for the global and local synchronization, and information flow capacity properties 
for the evolved networks of Materials and Methods, Subsection A Model for Brain Network Evo¬ 
lution Based on the Maximization of Information Flow Capacity. Panel (A): Parameter space for the 
synchronization (p ) 5 . Panel (B): Parameter space for the averaged upper bound for MIR, (mMIR) 5 , from the five 
realizations of a network of 60 neurons with six, equally sized, small-world clusters. Case A of high synchronization 
and low information flow capacity is denoted by A and case B of low synchronization and high information flow 
capacity by •. Panels (C) to (H) are plots for the local synchronization p Ci of the six communities of the C.elegans 
brain network. To be compared with panel (A). Here, g n is the chemical and gi the electrical coupling of Eqs. ([2]). 


with an initial small-world network configuration, by maximizing information flow capacity, the final network 
exhibits not only similar topological (structural network characteristics) but also similar functional or behavioral 
(synchronization and upper bound for MIR) properties as those found for the brain dynamical networks of the 
C.elegans and human subjects (see Materials and Methods, Subsection Spectral Similarity of C.elegans and Human 
Brain Networks with those of the Model for Brain Network Evolution). 

For completeness, in Figure S2, we present a similar analysis to the one of Fig. [3] based on Erdos-Renyi random 
networks (panels (A), (B)), scale-free (Barabasi-Albert) (panels (C), (D)) and star topologies perturbed by 20% for 
the clusters of the model for brain evolution of 60 neurons and 6 small-world clusters (panels (E), (F)) and found 
out that the evolution model fails to capture similar functional (local and global synchronization and, information 
flow capacity patterns in the parameter spaces) as the same model equipped with small-world topologies for its 
clusters (see Fig. |3|. Here, we are interested in studying and proposing a model for brain network evolution 
that is able to reproduce not only similar functional properties such as information flow capacity and, local and 
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global synchronization properties, but also importantly to reproduce similar structural properties for the finally 
evolved full brain network and of its clusters. Based on our results so far, we show next that the initial cluster 
configuration able to fulfil both requirements is the small-world cluster topology. 


Spectral Similarity of C.elegans and Human Brain Networks with those of the Model 
for Brain Network Evolution 

The discussion about the results of Fig. [4] in Materials and Methods suggest that the evolution process of a basic 
small-world clustered network is capable of generating an evolved one with similar structural properties with those 
for the C.elegans and human BDNs (see Materials and Methods, Subsection Structural Properties of the Model for 
Brain Network Evolution). The structural similarity to the human brain network is even more remarkable if the 
couplings of the network to be evolved are within the range that promotes high levels of information flow capacity 
and low neural synchronization, a prominent example of which is case B (for a definition and discussion about 
cases A, B 1 see Materials and Methods, Subsection Brain Network Evolution Promotes Global no Hebbian-like 
and, Local Hebbian-like and no Hebbian-like Evolution Learning). 

Here we study how close the normalized Laplacian spectral plots of the networks considered in Fig. BQ) are 
from those of evolved BDNs. Particularly, we examined the spectral similarity by comparing spectral plots and 
computing their average Euclidean distance following Ref. [35]. We refer the reader to Materials and Methods, 
Subsection Normalized Laplacian Spectra, for the details of our computations. We compared the normalized 
Laplacian spectral plot of the C.elegans and of the six human brain networks with the spectral plots of all finally 
evolved networks of the five realizations used to compute the averaged parameter spaces of the model for brain 
network evolution of Fig. §A), (B). We provide additional support by demonstrating in the last four panels of 
Fig- H results from a similar comparison between the C.elegans and humans with a double-sized version of 120 
neurons of the model for brain network evolution of the Materials and Methods, Subsection A Model for Brain 
Network Evolution Based on the Maximization of Information Flow Capacity. The normalized Laplacian spectra 
for the C.elegans and for the averaged over the six human BDNs are shown in Fig. BQ)- As it can be seen, both 
spectra share some common features: Both show a left-skewed distribution in which the largest eigenvalue is closer 
to one in agreement with results in Ref. 35 . Also, the distributions show peaks around one and the eigenvalues 


are scattered at the beginning of the spectra, suggesting similarities in their community structure, reminiscent of 
their small-worldness. 

Although the spectral plots of both cases of the model for brain network evolution shown in Fig. |4](R) do 
not exhibit all properties of the spectral plots of the C.elegans and humans of Fig. Bq), maybe due to the 
considerably smaller size of the former networks, they do exhibit interesting similarities: They are both left-skewed 
distributions with a peak around 1.3, being closer to 1 than to 2. We also observe low relative frequency eigenvalues 
at the beginning of both spectra. Both spectral properties suggest similarities in their community structure 35 


as well (i.e. their small-worldness). This is in accordance with the spectral plots of the averaged human and 
C.elegans brain networks of Fig. BQ)- We ar gue that this similarity comes from the small-worldness of the 
communities. The close relation between the normalized Laplacian spectra of Fig. [4] suggests the existence of 
common underlying structural properties of the neural networks of the C. elegans , the humans and the model for 
brain network evolution. 

We measured the similarity between the spectral plots of the C. elegans and the averaged human brain network 
with the averaged model for brain network evolution and, plot in panels (A), (B) of Fig. [5] their spectral distance 
D for different chemical and electrical coupling ranges. In this framework, the closer D is to zero, the closer 
structurally the compared networks are. In particular, panel (A) is the parameter space for the spectral distance 
between the C.elegans and the averaged model for brain network evolution and, panel (B) is a similar plot for the 
spectral distance between the averaged human brain network and the same model for brain network evolution. We 
pinpoint case A by a A and case B by •. 

Such results allow one to draw interesting relations between structure and function of the proposed model for 
brain network evolution with the structural properties of the brain networks of C.elegans and humans. Panels (A) 
and (B) of Fig. [5] already reveals that the smallest mean spectral distance happens for the pair of chemical and 
electrical couplings that gives rise to BDNs that present small amount of synchronization and high information 
flow capacity, in other words to cases such as B. In contrast, one of the largest mean spectral distances was found 
for case A that promotes high amount of neural synchronization and small information flow capacity in the brain 
network! 
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Figure 4: Structural properties and normalized Laplacian spectra of the brain networks considered 
in this study. Panels (A) to (D): Plot of the pdf of the normalized degrees ki (panel (A)), plot of the clustering 
coefficient CC(ki ) (panel (B)), plot of the average normalized degree k nn (ki) of the neighbors of nodes with 
normalized degree ki (panel (C)) and the network with its distinct clusters and communities given by different 
colors for the case A of high synchronization and low mMIR of the model for brain network evolution of 60 neurons 
and 6 clusters. Panels (E) to (H): Same as in panels (A) to (D) but for the case B of low synchronization and high 
mMIR of the same model. Panels (I) to (L): Same as in panels (A) to (D) but for the C.elegans brain network. 
Panels (M) to (P): Same as in panels (A) to (D) but for the human subject Ai brain network. In the plots of the 
second column, we show with blue dashed lines the exponential dependence of CC (ki) to ki to guide the eye, where 
ki is the normalized degree. It is defined as ki = ki/k max , where k t is the node degree and k max is the largest 
node degree in the network. For the first row, k max = 7, for the second k max = 8, for the third k max = 76 and for 
the fourth, k max = 87. Panel (Q): Normalized Laplacian spectra of the brain network of C.elegans (solid curve) 
and of the averaged over the six human subjects (dashed curve). Panel (R): Similarly for the brain network of case 
A of high synchronization and low information flow capacity (solid curve) and, for case B of low synchronization 
and high information flow capacity (dashed curve) of the model for brain network evolution. 
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Figure 5: Spectral distances between the network topology of the C.elegans , humans and evolved 
networks and, global synchronization and information flow capacity properties for an evolved BDN 
of 120 neurons. Panel (A): Parameter space for the spectral distance D between the C.elegans brain network and 
the averaged model for brain network evolution of 60 neurons and six small-world clusters and, panel (B) similarly 
for the spectral distance between the averaged brain network of the six human subjects and the same averaged 
network created by our brain network evolution. ▲ denotes case A and • case B, both explained in Materials and 
Methods, Subsection Brain Network Evolution Promotes Global no Hebbian-like and, Local Hebbian-like and no 
Hebbian-like Evolution Learning Processes. Both panels to be compared with Fig. (3^B). Panel (C): Parameter 
space for the synchronization p and panel (D) for the rnMIR of evolved networks using our brain network evolution 
process with 120 neurons. Panels (E), (F) are similar to (A), (B) for the spectral distance between the C.elegans , 
averaged brain network of the six human subjects and the model for brain network evolution of 120 neurons. 
Panels (E), (F) to be compared with panel (D). Here, g n is the chemical and gi the electrical coupling of Eqs. (§. 


10 



















Materials and Methods 


Preparation of Data for the Study of BDNs 


C.elegans Data C.elegans is the most commonly used model organism for neural network studies. It is a 1mm 
long soil worm with a simple nervous system that can be represented at first by 302 neurons and about 7000 
synapses [32) . Its nervous system is divided into two distinct and independent nervous systems: A large somatic 
nervous system with 282 neurons and a small pharyngeal nervous system with only 20 neurons. Since neurons 
CANL/R and VC06 do not make synapses with other neurons (i.e. they are isolated), they can be removed and a 
reduced connectivity brain matrix of 279 somatic neurons can be produced. We use in our study the connectome of 
the large somatic nervous system found in Ref. 36 that consists of 277 neurons. We decided to use the undirected 


version of this adjacency matrix as we are not concerned with the directionality of the information flow. We 
simulate the dynamics of each “neuron” by a single HR neuron system given in Eq. ([!]) and couple them by the 
corresponding adjacency matrix obtained from the brain connectivity of the C.elegans using Eqs. ([2]). 

We study the C.elegans nervous system in order to understand the human nervous system. The reason is that 
both human and C.elegans nervous systems consist of neurons and the communication or flow of information is 
passing through synapses that use neurotransmitters to perform brain activity. Many of these neurotransmitters 
are common in humans and C.elegans such as Glutamate, GABA, Acetylcholine and Dopamine. The genome of 
the C.elegans is almost 30 times smaller than that of humans but still, encodes almost 22000 proteins. Moreover, 
it is almost 35% similar to that of humans 


37 38 


Human Subjects Data The data for the analysis of the human connectome were based on Refs. 33 36 . The 


authors report on results based on the study of five different subjects (right handed males aged between 24 and 
32) coded as A, B, C, D and E, where the first one was examined twice giving connectomes A 1; A 2 with the 
second examination performed several days after the first one which yielded a highly consistent regional adjacency 
matrix A 2 to Ai. In our study, we make use of all six adjacency matrices (referring thereafter to six human 
subject BDNs) to prepare averaged quantities for global synchronization and upper bound for MIR. The diffusion 
spectrum imaging technique was then employed to retrieve high-resolution connection matrices for all subjects. 
The cortical regions of their brains were then further divided into 66 clear anatomical regions and these were 
individually subdivided into smaller regions of interest (of size 1.5cm 2 ) which finally resulted in 998 parts. These 
parts cover the entire cortices of both hemispheres but do not include subcortical nodes and connections 33 . Thus, 


the original adjacency matrices obtained for brain circuitry had 998 neural ensembles for all subjects. However, 
during the analysis of the adjacency matrices of these subjects, we found out that all these brain networks were 
disconnected and different numbers of isolated neural ensembles (nodes) were identified for each connectome. We 
decided to remove the isolated nodes as they were not connected to other neural ensembles of the connectome, 
and ended up with 994 for subject Ai, 987 for A 2 , 980 for B, 996 for C and D, and 992 for E. We also decided to 
use the undirected versions of the adjacency matrices for the same reason as for the C.elegans. We simulated each 
neural ensemble by a single HR neuron given in Eq. ([l]) and coupled them by the corresponding adjacency matrix 
obtained from the brain connectivity using Eqs. <§• 


Hindmarsh-Rose Neural Model for Brain Dynamics 

The complexity of the circuitry of the nervous system of the human brain is still a big challenge to be resolved as 
it contains about 86 billion neurons and thousands times more synapses 39 . A synapse is a junction between 


two neurons and it is a mean through which neurons communicate with each other. There are electrical and 
chemical synapses: An electrical synapse is a physical connection between two neurons which allows electrons to 
pass through neurons by a very small gap between nerve cells. Electrical synapses are bidirectional and of a local 
character, happening between neurons whose cells are close. They are believed to contribute to the regulation of 
synchronization in the brain network. In contrast, chemical synapses are special junctions through which the axon 
of the pre-synaptic neuron comes close to the post-synaptic cell membrane of another neuron or non-neural cell. 
In Ref. 40 , the authors report on the self-consistent gap junctions and chemical synapses in the connectome of 


the C.elegans. In our work, we use both kinds of synapses. 

Following Ref. 34 , we endow the nodes (i.e. neurons for the c.elegans and neural ensembles for the humans) 
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of the networks with Hindmarsh-Rose brain dynamics 41 


p = q-ap 3 + bp 2 -n + / ext , 
q = c - dp 2 - q, 
n = r[s(p - po) - n\, 


( 1 ) 


where p is the membrane potential, q is associated with the fast current, NaA or K + , and n with the slow current, 
for example Ca 2+ . The rest of the parameters are defined as a = 1, b = 3, c = 1, d = 5, s = 4, po = —1.6 
and I ex t = 3.25 for which the system exhibits a multi-scale chaotic behavior characterized as spike bursting. 
r modulates the slow dynamics of the system and was set to 0.005 so that each neuron or neural ensemble be 
chaotic. For these parameters, the HR model enables the spiking-bursting behavior of the membrane potential 
observed in experiments made with a single neuron in vitro. It is also a relatively simple model that provides a 
good qualitative description of the many different patterns empirically observed in neural activity. 

We couple the HR system to create an undirected BDN of N n neurons connected simultaneously by electrical 
(linear diffusive coupling) and chemical (nonlinear coupling) synapses: 


N n N n 

Pi = qi — ap 3 + bp 2 -rn + / ext - g n (pi - Vsyn) ^ ~ 91^2 

j =i j=i 

q i =c-dp 2 - qi, 
hi = r[s(pi -p 0 ) - Hi], 
i _ qiPi - Mi . Ar 

pf + qf 


( 2 ) 


In our study, fa is the instantaneous angular frequency of the i-tli neuron and fa is the phase defined by the fast 
variables ( Pi,qi ) of the i-th neuron. We consider H(pi) = pi and: 


S{Pj) = 


l _)_ e - A (Pj- e syn) : 


(3) 


with d S y n = —0.25, A = 10, and Vsyn = 2 to create excitatory BDNs. Equation ([3]) is a sigmoidal function that 
acts as a continuous mechanism for the activation and deactivation of the chemical synapses and, also allows for 
analytical calculations of the synchronous modes and synchronization manifolds of the coupled system of Eqs. 
§ [34]. In Eqs. §, g n is the coupling strength associated to the chemical synapses and gi to the electrical 
synapses. For the chosen parameters, we have \p. t \ < 2 and that {pi — Vsyn) is always negative for excitatory 
networks. If two neurons are connected under an excitatory synapse then, when the presynaptic neuron spikes, it 
induces the postsynaptic neuron to spike. We adopt only excitatory chemical synapses in our analysis. We use as 
initial conditions for each neuron i: pi = —1.30784489 + r?[, qi = —7.32183132 + i)\, rii = 3.35299859 + r/[ and 


fa = 0, where rfi is a uniformly distributed random number in [0,0.5] for all * = 1,..., N n , following Ref. 34 


These initial conditions place the trajectory quickly in the attractor of the dynamics and thus, there is less need to 
consider longer transients. 

Gij accounts for the way neurons are electrically (diffusively) coupled and it is a Laplacian matrix (i.e. 
Gij = Kij — A,j, where A is the binary adjacency matrix of the electrical connections and K is the degree 
identity matrix based on A), and so Gij = 0. By a binary adjacency matrix, we mean an adjacency matrix 

with entries either 0 (no connection) or 1 (connection). B,., is a binary adjacency matrix and describes how the 
neurons are chemically connected and therefore its diagonal elements are equal to 0, giving thus B ij = hi, 

where fa is the degree of the i -th neuron, i.e. it represents the number of chemical links that neuron i receives 
from all other j neurons in the network. A positive (i.e. 1) off-diagonal value in both matrices A , B in row i and 
column j means that neuron i perturbs neuron j with an intensity given by giGij (electrical diffusive coupling) 
or g n Bij (chemical excitatory coupling), respectively. Therefore, the binary adjacency matrices C of the BDNs 
considered in this work are given by: 

C = A + B. 


Numerical Simulations Details 

We numerically integrated Eqs. © in Fortran 90 using the Euler integration method (order one) with time step 
St = 0.01. We decided to do so to reduce the numerical complexity and CPU time of the required simulations 
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to feasible levels as a preliminary comparison of trajectories computed for the same parameters (i.e. St, initial 
conditions, etc.) with integration methods of order 2, 3 and 4, revealed similar results. 

We evolve the dynamics of the brain networks and calculate their two largest Lyapunov exponents Ai, A 2 for 
the estimation of the upper bound for MIR, I c . We use the well-known method of Refs. [42j[43] to compute the 
Lyapunov exponents needed for the estimation of the upper bound I c for MIR (see also Materials and Methods, 
Subsection Upper Bound for MIR). The numerical integration of the HR system of Eqs. ([2]) for the C.elegans 
and human BDNs was performed for the final integration time tf = 5000 and the computation of the different 
quantities such as the order parameter p of Materials and Methods, Subsection Synchronization Measures in 
BDNs and the Lyapunov exponents, starts after the transient time t t = 300 to make sure that orbits converged 
to the attractor of the dynamics. The same parameters for the model of brain network evolution of Materials 
and Methods, Subsection A Model for Brain Network Evolution Based on the Maximization of Information Flow 
Capacity were set to tf = 2500 and tt = 300 to reduce the numerical complexity and CPU time to feasible levels, 
retaining similar results. We have been careful to check that using exactly the same values as for the C.elegans 
and humans, the conclusions were practically the same. 


A Model for Brain Network Evolution Based on the Maximization of Information 
Flow Capacity 

In this work we propose an artificial evolution model for brain network connectivity that captures important 
structural and functional properties of the BDNs of the C.elegans and humans. Our idea is reminiscent of modular 
processors that are sufficiently isolated and dynamically differentiated to achieve independent computations, but 
also globally connected to be integrated in coherent functions 1 . 

The model is based on the consideration of the combined effect of chemical and electrical synapses between 
neurons based on a topology reminiscent of interconnected brain network communities found in the BDNs of 
C.elegans and humans. In this study, we consider chemical synapses solely for the communication of neurons 
of different clusters of the network (inter-cluster connections) and electrical synapses for the communication of 
neurons within each cluster (intra-cluster connections). This idea comes from the biological local and non-local 
nature of these connections. 

Particularly, we consider a starting network topology for brain network evolution, where N c clusters of 
electrically coupled neurons are connected in a closed ring as shown in Figure SI of the Supporting Information. 
We endow each cluster with a small-world topology 191 as this is what we found to be more plausible to happen 
on the brain networks of the C.elegans and humans (see Subsection Analysis of Networks and Communitie^] in 
Materials and Methods). We also use in our evolution model, for simplicity but this is not mandatory, clusters of 
the same number of neurons. We denote the total amount of neurons in the network by N n . Each small-world 
cluster is connected by an inter-cluster connection with its two nearest neighbour clusters by only chemical 
excitatory connections. In Figure SI of the Supporting Information, one can see such an example of a small-world 
network topology that comprises N c = 6 clusters and N n = 60 neurons, where the red links denote the chemical 
inter-cluster connections and black the electrical intra-cluster connections. For this model network, we subsequently 
compute the two largest Lyapunov exponents Ai, A 2 of the BDNs following Refs. 42 43 to estimate the upper 


bound I c = Ai — A 2 for the MIR of the network, i.e. the maximum amount of information per time unit that can 
be exchanged between the neurons of the basic (not yet evolved) network, aka its information flow capacity (for 
the details see Subsection Upper Bound for MIFlfl in Materials and Methods). 

We then evolve the starting network, such as the one in Figure SI, by adding new chemical excitatory 
inter-cluster connections to simulate the creation of new chemical synapses between neurons of different clusters. 
The electrical connections, topology and the values of the chemical and electrical coupling strengths are not 
modified during the evolution process. We adopt the following evolutionary rule to imitate brain plasticity 19 : If 


the newly added inter-cluster chemical connection leads to an increase of I c prior to the addition, the new synapse 
is retained. If, instead, it is found not to increase I c then it is deleted from the network and the random search for 
another one starts, being this procedure iterative. We choose the nodes of the different small-world clusters so 
that to simulate the addition of new inter-cluster connections in a random fashion (i.e. the candidate links are 
randomly chosen from a uniform distribution). 

The iterative procedure is repeated until the maximum number of possible pairs of neurons from different 
clusters is exhausted. We denote by mMIR the value of I c of the finally evolved BDN which is always bigger or 
equal than the I c of the starting BDN. For different values of the coupling strengths, mMIR can be achieved for 
different numbers of added interconnections. In all cases studied, we also compute the global synchronization 
measure p of the finally evolved BDN (for the details, see Materials and Methods, Subsection Synchronization 
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Measures in BDNs) to allow for direct comparisons with mMIR and for the identification of relations between 
synchronization and information flow capacity in the BDNs. 


Rescaling of Chemical and Electrical Couplings for Parameter Spaces of Networks 
with Different Eigenvalue Spectra 


Following Ref. 1341. 


a rough estimation on the range of chemical and electrical couplings based on those used for 

and 


the parameter space of another network that is capable of reproducing similar synchronous behaviors 28 


similar amounts of Kolmogorov-Sinai entropy [29] , can be computed as following: Suppose we have produced a 
parameter space such as those of Fig. [2] showing behaviors of the BDN of the C.elegans as a function of g n and 
gi. Let us denote the maximum electrical coupling as g p, the maximum chemical as g p, the smallest positive 
eigenvalue of the Laplacian matrix of the electrical connections as wp and the average degree of the chemical 
connections as d c . Suppose now we want to compute a similar parameter space for another BDN. Let us denote 
by <7/ max , and d max the corresponding values of the new parameter space. Then, for the new maximum 

couplings we have: 


9i 


d c 

^max 

, , C 

^max 


gf (chemical coupling), 


gf (electrical coupling). 


(4) 

(5) 


To arrive at Eqs. using the results of Ref. [34], we have assumed that a network with an average degree d 

for its chemical connections behaves similarly to a network with the same degree for its chemical connections. 

Equations Q and ([5]) provide a rough estimation on the maximum coupling strengths that can be used for the 
new parameter spaces. Table S3 presents w m , d, g" lax and g ; max for the different BDNs considered in our work. 
Based on these rough predictions, we then identified as best matching ranges, those depicted in the figures of the 
paper. Based on the maximum values of the parameter space ranges of Fig. §A) for the C.elegans (gf = 0.3 
and gf = 2), we get for the average humans g™ ax = 0.21 and g™ ax = 1.71, in good agreement with the maximum 
values of the ranges in Fig. [2jC) and (D). For the model for brain network evolution we get g™ ax = 1.14 to 2.28 
and g™ ax = 0.97 to 1.17 depending on the particular BDN, in accordance with the range of couplings used in Figs. 
[3] and, 2 in Figure S2. Similarly, for the large version of our model for brain network evolution ( N n = 120, N c = 6) 
we estimated g“ ax = 2.16 and g” lax = 2.9, consistent with the maximum coupling strengths used in the last four 
panels of Fig. [5] Consequently, our methodology allowed us to identify regions of synchronization p and upper 
bound for MIR, / c , for the different BDNs of this work with similar functional and structural properties. 


Analysis of Networks and Communities 

We initially identified the communities of the networks using the walktrap method [44] of the igraph software 
with six steps. The algorithm detects communities through a series of short random walks, with the idea that the 
vertices encountered on any given random walk are more likely to be within a community. The algorithm initially 
treats all nodes as communities of their own, then merges them into larger communities, and these into still larger, 
and so on. Essentially, it tries to find densely connected subgraphs (i.e. communities) in a graph via random 
walks. The idea is that short random walks tend to stay in the same community. Following this procedure we 
have been able to identify 6 communities in the C.elegans BDN, 10 in human subject Ai, 5 in A 2 , 9 in B, 6 in C, 
10 in D and finally, 7 in E. 

After this step, we computed various statistical quantities such as the global clustering coefficient, the average 
of local clustering coefficients, the mean shortest path, the degree pdf of the network and the small-worldness 
measure. The latter property is characterized by a relatively short minimum path length on average between all 
pairs of nodes in the network, together with a high clustering coefficient. 

Even though small-worldness captures important aspects of complex networks at the local and global scale of 
the structure, it does not provide information about the intermediate scale. Properties of the intermediate scale 


can be more completely described by the community structure or modularity of the network 45 . The modules of 


a complex network, also-called communities, are subsets of nodes that are densely connected to other nodes in 
the same module but sparsely connected to nodes belonging to other communities. Since nodes within the same 
module are densely intra-connected, the number of triangles in a modular network is larger than in a random graph 
of the same size and degree distribution, while the existence of a few links between nodes in different modules 
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plays the role of topological shortcuts in the small-world topology. Systems characterized by this property tend to 
be small-world networks, with high clustering coefficient and short path length with respect to random networks. 

To infer the small-worldness of a network or community, we first compute the mean local clustering coefficient 
C of the network or community and the mean of the local clustering coefficients of one hundred randomly created 
networks (CV)ioo of the same degree pdf with the studied network and also, the mean shortest path of the studied 
network L and the average of the mean shortest paths of the same one hundred random networks (T r )ioo- Watts 
and Strogatz [9| measured that many real-world networks have an average shortest path length comparable to those 
of a random network (L ~ L r ), and at the same time a clustering coefficient significantly higher than expected by 
random chance (C C r ). Then, they proposed a novel graph model, currently named the Watts-Strogatz model, 
with a small average shortest path length L , and a large clustering coefficient C. We adopt this as a working 
definition of a small-world network or community. Therefore, small-world networks are in between the limit cases 
of regular graphs with large L and C and random networks with small L and C. To quantify small-worldness we 
use the ratios 46 : 

- L - c tn 

^ {L r ) 100 ’ ^ (Cr) 100 ’ 

in such a way that, for a small-world network or community, we compute: 


cr = — > 1, 


(7) 


being the small-worldness measure. The higher is a from unity for a given network or community, the better it 
displays the small-world property. 

For completeness, we note that for the human subject D, the walk-trap community analysis with step equal to 
six detected eleven communities for which the last one comprised only one neuron. Hence, in all computations, we 
disregarded this community as a trivial case. We present the results of the above analysis in Table S4. We have 
performed all structural analyses of this paper using the igraph software. The evolution of the basic network of 
Figure SI, under the principle of the maximization of the information flow capacity, is able to capture behaviors of 
real brain connectivity networks such as those for the communities of the C.elegans BDN, being their small-world 
structure a prominent reason. We found out that in all networks studied, the small-world measure cr gets values 
much higher than unity (see Table S4), clearly indicating that they all display the small-world property, though in 
different degrees. 


Brain Network Evolution Promotes Global no Hebbian-like and, Local Hebbian-like 
and no Hebbian-like Evolution Learning Processes 

Here, based on the extraordinary functional similarities found so far, we focus on two characteristic behaviors of 
the model for brain network evolution defined in Materials and Methods, Subsection A Model for Brain Network 
Evolution Based on the Maximization of Information Flow Capacity, to reveal important modular behaviors and 
underlying learning evolution processes: We associate the first one to the case of high global synchronization 
and low information flow capacity (which we call case A) and the other one to the opposite situation of low 
synchronization and high information flow capacity (case B). For illustration purposes, we select from the proposed 
model for the brain network evolution of 60 neurons and six clusters of Materials and Methods, Subsection A 
Model for Brain Network Evolution Based on the Maximization of Information Flow Capacity, two finally evolved 
BDNs with the following coupling strengths: For case A the pair g n = 0.2, gi = 1.8 from the fifth realization 
(indicated by ▲ in Figs. [3j[5]) and for case B the pair g n = 0.9, gi = 1.5 from the fourth realization (indicated by • 
in Figs. Hi- We have checked that the conclusions are valid independently of the realization and other similar 
pairs of coupling strengths. Also, our results reported here are independent on the initial small-world cluster 
configuration and initial conditions. 

In panels (A), (B) of Fig. [6]we demonstrate the relation between global synchronization p and mMIR for these 
cases. Case B of moderately low global synchronization (dashed black curve in panel (A)) and high information 
flow capacity (dashed black curve in panel (B)) is characterized by a larger number of added interconnections (i.e. 
30) present in the finally evolved network with respect to case A of only 12 (corresponding to the solid black curves 
in panels (A), (B))! In both cases, p of Eq. (|8j) for global synchronous behavior, has the tendency to decrease with 
different slopes (denoted by 9 in Fig. |6](A)) during brain network evolution as is evident by the fitting to the data 
in dashed blue lines. We attribute this behavior to an evolutionary brain network behavior, where global neural 
synchronization levels decrease during brain network evolution, reminiscent of a possible underlying global no 
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Figure 6: Brain network evolution promotes Hebbian-like and no Hebbian-like processes and, mod¬ 
ular organization in the brain dynamical networks. Panels (A), (B): Global synchronization p and mMIR 
respectively for cases A and B. 9s are the slopes of the two dashed blue lines which are fitted to the black curves 
to demonstrate the decrease of the global neural synchrony during brain network evolution. Panels (C), (D) show 
how the average pair-wise synchronization ( pij) Cl of the clusters Ci, l = 1,... ,6 changes during evolution, for cases 
A, B respectively. Panels (E), (F) show the pair-wise neural synchronization level p,j of the finally evolved BDNs 
for the same cases. The horizontal axes (Number of links) correspond to the added links during brain network 
evolution that lead to the increase of the information flow capacity at each step. The caption of panel (D) for the 
black curves is the same as in panel (C). These results are for the studied model of brain network evolution with 
60 neurons and six small-world clusters. 


Hebbian-like evolution process that promotes a decrease in global synchronization levels as new connections are 
added in the network. 

We regard this as an important global property of the model for brain network evolution that needs to be 
further clarified, as we need to account for what happens on the local, cluster level as well. We thus use the 
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following procedure to infer about the underlying learning rules for the clusters: During each step of the evolution 
process of the model of 60 neurons and six small-world clusters (which is a particular BDN), we compute ( pij) ci 
based on Eq. © in Materials and Methods, to account for the average pair-wise synchronization of cluster 
Ci, l = 1,... ,6. At the end of the time evolution, we have six such values for the clusters, lying in the interval [0,1]. 
We then record these values if the newly added chemical interconnection leads to an increase of the information 
flow capacity as depicted by the corresponding I c and disregard them if not. At the end of the process, we 
result with a relationship between ( Pij) C i and the added interconnections that maximizes J c , for all clusters (see 
panels (C), (D) of Fig. [6]). Doing so, we can account for the underlying cluster learning processes. In particular, 
the terminology “Hebbian-like” is employed to represent a learning rule that not only involves quantities for 
synchronous events defined between pairs of neurons (see Eq. ©)> but also a direct relation between synapse 
strength and synchronization increase. The “no Hebbian-like” terminology refers to a learning rule that involves 
not only a global measure of synchronous behavior (see Eq. but also refers to a phenomenon where the 

addition of synapses is accompanied by a decrease in the synchronization levels. 

We present these results in Fig. [6](C), (D). Following our considerations, we already know that globally, both 
finally evolved BDNs are following a no Hebbian-like learning rule. However, panels (C), (D) reveal a substantial 
difference for the synchronization behavior for pairs of neurons in the clusters in the two cases that allows us to 
assign different kinds of learning rules to them. For case B in panel (D), the slopes 9 CI of the fitted lines (in blue) 
to the data ( pij ) ci , l = 1,... ,6 (in black) show that the third and fifth cluster have the tendency to increase their 
internal synchronization level as the BDN evolves. In other words, neurons belonging to these two clusters follow a 
Hebbian-like learning process that comes in contrast to the no Hebbian-like learning behavior of neurons belonging 
to the other clusters as they show the opposite trend during brain network evolution, exhibiting negative slopes 
for the synchronization. These findings are in agreement with the conclusions drawn from panel (F) in which 
neurons belonging to the third and fifth cluster are seen to be more synchronized with respect to the others (in 
yellow and red), a situation that promotes the modular organization in the brain with different internal levels of 
synchronization. On the other hand, case A demonstrates a completely different situation in which the level of 
global synchronization decreases and cluster synchronization is very strong between neurons in all clusters during 
the whole brain network evolution, a behavior that results in a network whose neurons are unable to exchange 
but very small amounts of information. This is a case where the brain network can not transmit information by 
increasing its synaptic efficacy and thus, it is unable to learn new information! 

Panels (E), (F) of Fig. [6] show the pair-wise synchronization level of the neurons denoted as pij (see Eq. © in 
Materials and Methods, Subsection Synchronization Measures in BDNs) of the finally evolved BDNs for cases A 
and B , resulting in networks of 12 and 30 chemical inter-links, respectively. Panel (F) of case B for moderately low 
global synchronization and high information flow capacity reveals a clusterized synchronization behavior. The 
third and fifth small-world clusters show much higher levels of internal synchrony (yellow and red points), i.e. 
synchronization of pairs of neurons in the same small-world community, with respect to the blue or dark blue 
points of low neural pair synchronization in other small-world clusters. The same also happens for a small number 
of pairs of neurons belonging to different clusters (off-diagonal points). One can notice that given the cluster that 
has the same internal level of synchronization, there can always be found two subnetworks, each belonging to one 
cluster, that have an equal amount of synchronization. This means that case B corresponds to a BDN that has 
clustered multilayer synchronization, where different clusters become functionally connected with different common 
behaviors. A situation reminiscent of findings in neuroscience that demonstrate the modular organization of the 
brain in which modular processors are sufficiently isolated and dynamically differentiated to achieve independent 
computations, but also globally connected to be integrated in coherent functions 1 . Our results for case B (see 
Fig. [6|F)) come in contrast to those of A of high global synchronization and low information flow capacity shown 
in Fig. [6](E) which demonstrates that all clusters and pairs of neurons attain almost the same state of almost 
complete synchronization, revealing a highly synchronized brain dynamical network that is not able however to 
exchange but only very small amounts of information between its different parts! 

The no Hebbian-like mechanism is effectively similar to the unlearning anti-Hebbian mechanism of Crick and 
Mitchison 27 that proposes the elimination of unnecessary connections to prevent overload and to render the 


network more efficient. Both mechanisms lead however to more efficient networks, being in our case the evolved 
networks able to maximize their information flow capacity. 

Typically, low synchronization implies Ai ss A 2 leading to I c ss 0. In critical points however, I c tends to be large 
as long as most of the oscillation modes are stable resulting in a situation where Ai > 0 and, A 2 ~ 0 and positive. 
Thus, maximum I c at low synchronization corresponds to maximum I c near the critical point of the dynamics 
where Ai > 0 and, A 2 ~ 0 and positive. Self-critical phenomena happen when the network has marginal Lyapunov 
exponents, in other words at the critical point that results in typically large I c . In this context, our results are in 
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agreement with the work in Ref. 47 , supporting our findings for the existence of Hebbian and no Hebbian-like 


learning mechanisms in the level of the communities (local synchronization) at self-criticality, which is responsible 
for the maximization of the information flow capacity of the evolved BDNs (low global synchronization), such as 
in case B (for a similar result see Ref. [48]). 

We do not study our BDNs under different stimuli. Our main hypothesis is that the final topology and behavior 
of an evolved BDN that maximizes MIR between its neurons is similar to real brain networks, such as those from 
the C.elegans and human subjects. Since our results for the maximization of the information flow capacity of the 
evolved BDNs happen when self-critical phenomena emerge, they are in agreement with the results reported in 
Refs. 
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Structural Properties of the Model for Brain Network Evolution 

We study here the structural properties of the evolved BDNs and compare with those of the C.elegans and human 
BDNs. As we have already demonstrated, they share common functional properties. 

We present the results of this study in Fig. [4] Particularly, panel (A) shows the degree probability distribution 
function (pdf(fcj)), panel (B) the clustering coefficient CC (ki) as a function of the normalized degree ki , panel (C) 
the average degree k nn (ki) of the neighbors of nodes with degree ki and panel (D) the network with its distinct 
communities depicted by different color-shaded neighborhoods, for case A. In this context, ki is the normalized 
with respect to the maximum, degree. Panels (E) to (H) show similar plots for case B. Panels (I) to (L) are similar 
plots for the C.elegans brain network and, panels (M) to (P) for human subject Ai. The plots in the second 
column show the correlation between different normalized degrees of the network whereas those of the third the 
tendency of the nodes of a certain normalized degree k to link with other nodes of a given degree. From the plots 
of the second column we observe that high degree nodes have the tendency to link with low degree nodes following 
an exponential dependence (with different exponents) implying disassortative mixing by degree. Our results from 
the second column of Fig. [4] suggest that evolving BDNs based on the maximization of the upper bound for MIR 
(cases A and B) gives rise to disassortative mixing by degree meaning that high degree nodes are preferentially 
connected to other low degree nodes and low to high degree nodes. 

It is worth noting that the structural properties of all human subjects are similar. Figures |4](I) to (O) show 
that human subject Ai has similar structural properties to the C.elegans structure. Cases A and B seem also 
to present strong structural similarities, for the quantities considered in Fig. [4] despite the profoundly different 
functional behaviors. We found out that case B of moderately low global synchronization and high information 
flow capacity is characterized by a big number of added interconnections (i.e. 30) present in the final network 
whereas case A by only 12! We have checked that this relationship between a large (small) number of added 
chemical inter-links with low (high) synchronization and large (low) information flow capacity happens for all 
similar cases of the parameter spaces. It is also noteworthy that in case B (see Fig. ]4|H)), the number of clusters 
of the finally evolved network is reduced by one (since we start with six and end up with five) and is in contrast 
with what is happening in case AI Depending on the couplings which give rise to different functional behaviors, 
the brain network evolution model is capable of merging different communities, i.e. of restructuring the initial 
network configuration. Case B is a more globally connected network. 

The previous observation can be quantified in terms of the modularity of the network, i.e. by the strength of 
the division of the network into modules (groups, clusters or communities). Networks with high modularity have 
dense connections between the nodes within modules and sparse connections between nodes of different modules, 
being more clustered. We have used the igraph software for these computations. By applying this idea here, we 
find that the modularity of the final BDN of case B is 0.596, very close to the average modularity of the six human 
subjects 0.588 ± 0.023. In contrast, the modularity of case A is 0.702. For the sake of completeness, we also report 
the modularity of the C.elegans topology which is 0.375, the smaller of all cases we considered. The last result 
shows a network with sparser connections between the nodes within the modules and denser between nodes of 
different modules! For the C.elegans , the natural evolution process led to a smaller clusterization of its brain 
network. 

The results of Fig. [4] suggest that the evolution process of a basic small-world clustered network is capable 
of generating evolved ones with similar structural properties to those for the C.elegans and human BDNs. The 
structural similarity to the human brain is even more remarkable if the couplings of the network to be evolved are 
within the range that promotes high levels of information flow capacity such as in case B. 
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Synchronization Measures in BDNs 


Synchronous activity has been observed in neural systems and reported to be associated not only with pathological 
brain states 


51 but also with various cognitive functions 52 


It has been found that burst synchronization of 

noise 


neural systems may be strongly influenced by many factors, such as coupling strengths and types 
and the existence of clusters in neural networks 


14 


53 


54 


In this paper we use the order parameter p to account for the synchronization level of the neural activity of the 


studied BDNs and of their communities 55 . It is originated from the theory of measures of dynamical coherence 
of a population of N n oscillators of the Kuramoto type 56 and, can be computed by a complex number z(t) 
defined as: 


Z {t) = P (ty*w = 


N n 

y] (*) 

3 = 1 


( 8 ) 


By taking the modulus p(t) of z(t), one can measure the phase coherence of the population of the N n neurons of 
the BDN, and by $(t) to measure the average phase of the population of oscillators. In this context, <fii is the 
phase variable of the *-th neuron of the HR system ([2]) given by its fourth equation. Actually, one averages over 
time p(t) to obtain the order parameter p = ( p{t )) t , the tendency of p in time. A value of p = 1 corresponds to 
complete synchronization of the oscillators, whereas p = 0 to complete desynchronization. 

We use Eq. ([8]), adapted accordingly, wherever in the paper we need to compute the synchronization level 
of BDNs or clusters. In particular, in the case of BDNs, N n is the number of neurons of the BDN and j runs 
through all N n neurons of that network whereas in the case of clusters, N n represents the number of neurons of 
the particular cluster and j refers to the particular neurons which are members of this cluster. 

We also compute and plot in Fig. [6^C), (D), for cases A and B respectively of the model for brain network 
evolution of 60 neurons and six small-world clusters, the pair-wise neural synchronization by looking at the 
synchronization patterns between all pairs of neurons i, j of the network as: 


Pij — 


lim 
At—too 


s-i nr+At 

/ gi (*)](/£ 

A t L 


(9) 


where C % j is the adjacency matrix of the brain network and <j>i is the phase variable of the i-th neuron of system 
([2|. p^ are bounded in the interval [0,1], being p^ = 1 when neurons i, j are fully synchronized and 0 when they 
are dynamically uncorrelated. To correctly compute p ^, we take the averaging time large enough in order to 
obtain good measurements of the coherence degree of each pair. We similarly compute and plot in Fig. §C), 
(D) for the same cases, the averaged quantity ( pij)cu / = 1,. - -, 6 over the six clusters, where i,j run through the 
neurons of each cluster. 


Upper Bound for MIR 


After Shannon’s pioneering work 57 on information, it became clear that it is a very useful and important concept 


as it can measure the amount of uncertainty an observer has about a random event and thus provides a measure 
of how unpredictable it is. Another related concept to the Shannon entropy that can characterize random complex 


systems is the MI 57 which is a measure of how much uncertainty one has about a state variable after observing 


another state variable. 

In Ref. 25 , the authors have derived an upper bound for the MIR between two nodes or two groups of nodes 
of a complex dynamical network that depends on the two largest Lyapunov exponents l\, 


network formed by these nodes. In particular, they have shown that: 

MIR < t I r = / j — l“2 , ^1 A 


I 2 of the subspace of the 


( 10 ) 


where Zi, I 2 are the two finite time and size Lyapunov exponents calculated in the bi-dimensional observation 
space of the two considered nodes 25 58 , which typically should approach the two largest Lyapunov exponents 
Ai, A 2 of the dynamical network if it is connected and the time considered to calculate l\, I 2 is sufficiently small. 
In our study, the upper bound I c for the MIR between any two nodes of the BDNs is effectively estimated by 
I c = Ai — A 2 (i.e. l\ = Ai and I 2 = X 2 ) and will stand for the upper bound for the information transferred per 
time unit between any two nodes of the BDN (i.e. between the neurons), what represents the information flow 
capacity of the BDN. The phase spaces of the dynamical systems associated to the neural networks we study 
here are excessively highly multi-dimensional and thus, estimating an upper bound for the MIR using Ai and A 2 
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calculated by the methods of Refs. [42] 43 instead of the MIR itself between all pairs of nodes, reduces enormously 
the computational complexity of the numerical calculations of this work. Besides, parameter changes that causes 
positive or negative changes in the MIR are reflected in the upper bound with the same proportion 


25 


Normalized Laplacian Spectra 


It is well-known that similarities between the structure of networks can be used for their classification 59 


The architecture of brain networks that describe the organization of maps of connections between neurons and 
brain elements at a systems level can be achieved by examining the eigenvalue spectrum of the normalized 
Laplacian of the connectome 35 60 . In our study, the connectome is given by the adjacency matrix of the brain 


network and thus we compute the eigenvalues of the normalized Laplacian based on this matrix. The eigenvalues 
Vi, i = 1,..., N n of the normalized Laplacian matrix L are in [0,2] which helps to compare networks of different 
sizes, and is defined as: 

if i = j, 

if i,j are connected, (11) 

otherwise, 


1 


Lij — 


0 


where i,j represent any two nodes of the brain network, L l3 the link between nodes i, j and ki the degree of node 
i. Therefore, the Laplacian spectrum of the network is given by the set of all eigenvalues of L, namely by its 
eigenspectrum. 

Spectral Plot The spectral plots were obtained from a smoothed eigenvalue distribution T(a;) that consists of 
eigenvalue frequencies convolved with a Gaussian kernel 
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where N n is the number of eigenvalues of L and a a smoothing function. In our study, we used a = 0.015 to 
smooth out appropriately the spectral plots. For these plots, a discrete smoothed spectrum was used in which T 
had steps of 0.001 and the distribution was normalized in such a way that the total eigenvalue frequency is unity. 


59 , and is defined as the average Euclidean 


Spectral Graph Distance The similarity distance between spectral plots was quantified using a spectral 
distance measure, based on the distance measure introduced in Ref. 
distance between two spectral plots Ti and r 2 : 

1 


m,r 2 ) = —^minfvM^WTHF) 

E“ in (V[ri«-r 2 o -)] 2 + [i-j] 2 ), 

j—0 1 ' ' 


k + 1 
1 

k + 1 


(13) 


where T(?’) is the discrete, normalized and smoothed eigenvalue distribution of Eq. (12) and the number of intervals 
k was set to 2000. The distance function (13) depends on the scaling of the axes, meaning that different scales 


result in different distances. Therefore, it is not an invariant distance measure between two networks but only 
serves as a tool to underpin the visual results in a quantitative manner. 

To generalize the important relation between function and structure in BDNs, we plot in panels (C) to (F) 
in Fig. [5] the results of a similar study for one realization of a double-sized model for brain network evolution 
based on 120 HR coupled neurons arranged in six small-world clusters. Panels (C) and (D) show the parameter 
spaces for the global synchronization p and information flow capacity mMIR of the big BDN which reproduce 
well the behavior of the small version of the model for brain network evolution of Fig. §A), (B) and allow for 
similar functional conclusions to be drawn. The relation between structural and functional properties for the 
double-sized model can be seen in the last two panels where we present the parameter spaces for the spectral 
distances when compared with the C.elegans (panel (E)) and averaged humans (panel (F)). Again, couplings that 
promote low global neural synchronization and high information flow capacity give rise to the biggest possible 
spectral similarity of the double-sized brain network evolution model with the C.elegans and human brain networks, 
the same conclusion drawn for the small version of the same model! We therefore verify our hypothesis that finally 
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evolved BDNs with neurons that can potentially exchange the highest levels of information are the BDNs with 
topologies closer to the brain networks of the C.elegans and humans. 


Discussion 

In this paper we propose a working hypothesis, and provide evidence, that neural networks that evolve based on 
the principle of the maximization of their internal information flow capabilities produce networks whose functional 
behavior and topology are similar to those features observed in dynamical neural networks whose topology is 
provided by the C.elegans and humans. 

Our hypothesis goes along the lines of the infomax theory that proposes that the brain evolves by maximizing 
the mutual information between external stimuli and its response. When maximizing the internal information 
flow capacity, we are creating a network capable of processing information about external stimuli for which its 
information content is smaller than the information flow capacity of the evolved network. Notably, the brain 
evolves by the action of input signals. Our working hypothesis simplifies enormously the complexity of the involved 
calculations and, allow us to understand function and behavior in the brain, without the need to externally perturb 
non-autonomous neural networks. 

We have been able to show that our evolved brain networks present similar synchronization and information 
flow capacity behaviors with the ones found for the simulated dynamical networks for the structure of the C. elegans 
and human brain. Moreover, we have shown that BDNs evolved with coupling strengths that maximize the 
information flow capacity are the ones that have the smallest spectral graph distance from the BDNs of the 
C.elegans and humans, and that, during the growing process, their MIR increase is related to moderately low 
amounts of global neural synchronization. Actually, the global neural synchronization levels decrease during brain 
network evolution, revealing an underlying global no Hebbian-like evolution process (where synapse strength leads 
to global decay of synchronization) driven by a mix of local no Hebbian-like learning rules for neurons in some 
clusters and by Hebbian-like learning rules in neurons belonging to other clusters where synapse strength leads to 
cluster synchronization. In this context, the no Hebbian-like mechanism is effectively similar to the unlearning 
anti-Hebbian mechanism as both lead to more efficient networks, in the sense that in our case the evolved networks 
are able to maximize their information flow capacity. 

We note that if other models than Hindmarsh-Rose will be used, such as the Morris-Lecar, Izhikevich, or spiking 
map-based neural models, then the parameter regions that maximize information flow capacity and minimize 
synchronization (or vice versa) will be different, however we expect that this would not change the main results and 
conclusions of this work in the sense that brain network evolution based on the maximization of information flow 
capacity will lead to similar topologies, behaviors and relations for the evolved networks. For the human subjects, 
the graphs represent functionally connected brain regions. The Wilson-Cowan model could be appropriate to 
model the human brain, whereas it will not be suitable to model the C.elegans brain. As we have done here, using 
neurons to represent nodes in the human connectome do not reproduce the real dynamics of the brain but gives us 
a mean to compare results with the C.elegans and the evolved networks. 

In relation to our work, the maximization of the information flow capacity for low synchronization corresponds 
to the critical point of the dynamics in which self-critical phenomena occur when the second or larger Lyapunov 
exponents of the BDNs are marginally positive. 

Finally, our results support further the hypothesis made in Ref. [2j that maximization of the information flow 
capacity can serve as a principle for the development of heterogeneous structures in brain dynamical networks, 
such as the neocortex of mammalian brains. 
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Figure SI: An example (first of the five realizations) of a starting small-world network topology for 
the proposed brain network evolution model. It comprises N n = 60 neurons arranged in N c = 6 small-world 
clusters. The red ring consists of chemical excitatory connections that link all clusters of the network. Within 
each cluster, depicted by a differently color-shaded area of intra-connected neurons, we consider solely electrical 
connections denoted by black edges. 


1 











S2 


A 

2 
1.5 
c5T 1 
0.5 
0 



1 

0.8 

0.6 

0.4 

0.2 

0 


0 0.2 0.4 0.6 0.8 1 

9n 



2 

1.5 

D) 1 

0.5 

0 



1 

0.8 

0.6 

0.4 

0.2 

0 


0 0.2 0.4 0.6 0.8 1 

9n 



2 
1.5 
d> 1 
0.5 
0 



1 

0.8 

0.6 

0.4 

0.2 

0 


0 0.2 0.4 0.6 0.8 1 



Figure S2: Global synchronization and information flow capacity properties for the evolved BDNs 
for different initial community structures. (A): Parameter space for the synchronization p and (B): for 
mMIR for six totally random (Erdos-Renyi) clusters. (C) and (D) are similar plots but for six scale-free (Barabasi- 
Albert) clusters. (E) and (F) are for six star clusters perturbed by 20%. All plots are from the model for brain 
network evolution of 60 neurons and six clusters. Here, g n is the chemical and gi the electrical coupling of Eqs. (2) 
of the paper. 
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Table S3: Smallest positive eigenvalue u> m of the Laplacian matrix G of the electrical connections 
and average degree of chemical connections d of the adjacency matrix B for the BDNs considered 
in this work. These values were used to provide rough estimates for the extend of the couplings of the parameter 
spaces, g,“ ax and g“ ax , based on those of the zoomed-in parameter space of the C.elegans of Fig. 2 of the main 
manuscript (first row of the Table). 
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Case A ( N n = 60, N c = 6) 
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Case B (N n = 60, N c = 6) 
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N n = 120, N c = 6 
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Average over the six human subjects 
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Table S4: Small-worldness measure <x for the C.elegans and human brain networks and for their 
communities. _ 
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